1 Demographics of H. schistosus across time

1.1 Sample size

snakes%>%
  group_by(Species)%>%
  count(name = "N")
Species N
Hydrophis curtus 282
Hydrophis schistosus 967
snakes%>%
  filter(Gear.Type != "")%>%
  select(Date, Gear.Type, No..of.Hauls, Average.Haul.Duration..Hours., Tow.duration.hours)%>%
  distinct()%>%
  mutate(Tow.duration.hours = case_when(!is.na(Tow.duration.hours) ~ Tow.duration.hours,
                                        is.na(Tow.duration.hours) ~ No..of.Hauls*Average.Haul.Duration..Hours.))%>%
  group_by(Gear.Type)%>%
  summarise(n = n(),
            haul.hours = sum(Tow.duration.hours, na.rm = T))
Gear.Type n haul.hours
GillNet 340 535.9633
Rampan 46 190.4500
Trawler 76 285.1600

1.2 Age structure

age.yr <- snakes%>%
  group_by(Species, Year)%>%
  summarise(N = n(),
            mean = mean(Snout.to.Vent..cm., na.rm = T))


maturtity <- snakes%>%
  group_by(Species)%>%
  count()%>%
  mutate(juv = 35,
         adult = ifelse(Species == "Hydrophis curtus", 54, 65))
# plotting distribution of SVL across years

snakes%>%
  filter(Snout.to.Vent..cm. > 20)%>%
  ggplot(aes(Year, Snout.to.Vent..cm.))+
  geom_violin(fill = "grey")+
  geom_boxplot(width = 0.1)+
  geom_hline(data = maturtity, aes(yintercept = adult), linetype = "dashed")+
  geom_hline(data = maturtity, aes(yintercept = juv), linetype = "dotted")+
  #stat_summary(fun.data = "mean_sdl", geom = "pointrange", size = 1)+
  geom_label(data = age.yr, aes(Year, 10, label = N))+
  facet_wrap(~Species, ncol = 1, scales = "free_y")+
  labs(y = "Snout to vent length (cm)")+
  theme(strip.text = element_text(face = "italic"))

ggsave(last_plot(), filename = "./Figures/figure1.tiff", height = 6, width = 8)
pop.yr <- snakes%>%
  filter(!is.na(Class))%>% # fix 2016 svl
  group_by(Species, Year)%>%
  count(Class)%>%
  complete(Class, fill = list(n = 0))%>%
  mutate(N = sum(n))%>%
  group_by(Year, Species, Class)%>%
  mutate(prop.age = n/N)

pop.yr%>%
  group_by(Species, Class)%>%
  skimr::skim(prop.age)%>%
  skimr::yank("numeric")%>%
  select(-hist)

Variable type: numeric

skim_variable Species Class n_missing complete_rate mean sd p0 p25 p50 p75 p100
prop.age Hydrophis curtus Adult 0 1 0.43 0.31 0.17 0.28 0.33 0.49 0.88
prop.age Hydrophis curtus Juvenile 0 1 0.57 0.31 0.12 0.50 0.66 0.72 0.83
prop.age Hydrophis curtus Neonate 0 1 0.02 NA 0.02 0.02 0.02 0.02 0.02
prop.age Hydrophis schistosus Adult 0 1 0.88 0.09 0.79 0.84 0.87 0.91 1.00
prop.age Hydrophis schistosus Juvenile 0 1 0.14 0.07 0.08 0.10 0.12 0.17 0.21
prop.age Hydrophis schistosus Neonate 0 1 0.03 0.02 0.01 0.02 0.03 0.03 0.04
pop.yr%>%
  filter(Year != 2019)%>%
  group_by(Species, Class)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)),
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(Class:p.value)
Species Class estimate1 estimate2 estimate3 statistic p.value
Hydrophis curtus Adult 0.8750000 0.3125000 0.3551402 9.2546983 0.0097807
Hydrophis curtus Juvenile 0.1250000 0.6875000 0.6261682 8.8509715 0.0119684
Hydrophis curtus Neonate NA NA NA 97.2336449 0.0000000
Hydrophis schistosus Adult 0.8846154 0.7894737 0.8613445 5.9949721 0.0499124
Hydrophis schistosus Juvenile 0.0769231 0.2105263 0.1239496 9.7219255 0.0077430
Hydrophis schistosus Neonate 0.0384615 0.0147059 NA 0.0189777 0.8904305

Age structure does not change significantly over a four year period from 2016 to 2019.

library(car)

snakes%>%
  filter(Year != "2019")%>%
  group_by(Species)%>%
  select(Year, Snout.to.Vent..cm.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
          sumr = map(mod, broom::tidy),
          stat = map(mod, car::Anova))%>%
  select(stat)%>%
  unnest()
Species Sum Sq Df F value Pr(>F)
Hydrophis schistosus 9866.294 2 21.86203 0.00e+00
Hydrophis schistosus 168785.548 748 NA NA
Hydrophis curtus 2641.334 2 10.07170 7.05e-05
Hydrophis curtus 24258.408 185 NA NA

GLM shows significance but how do you interpret these results?

# plotting distribution of SVL across months

month.svl <- HS%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month)%>%
  group_by(Month)%>%
  summarise(mean.SVL = mean(Snout.to.Vent..cm., na.rm = T))

births <- data.frame(Species = "Hydrophis schistosus", Month = "April")
  
HS%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(nesting(Species), Month)%>%
  droplevels()%>%
  ggplot(aes(Month, Snout.to.Vent..cm.))+
  geom_violin(fill = "light grey")+
  #geom_point(data = month.svl, aes(x = Month, y = mean.SVL), size = 3)+
  geom_boxplot(width = 0.1)+
  geom_segment(data = births, aes(x = Month, xend = Month, y = 0 , yend = 10), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(data = births, aes(x = Month, y = 20, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 80, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Snout to vent length (cm)")

Proportion of neonates increases in March to May. Birth and neonates were observed around the same time.

  • GLM testing change in proportion of SVL < 40 over months
  • Other ways to test age structure?
# proportion of jeuveniles in each month

snakes%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  group_by(Species, Month)%>%
  summarise(prop.neonate = sum(Snout.to.Vent..cm. < 40)/n())%>%
  spread(Month, prop.neonate)
Species January February March April May October November December
Hydrophis curtus 0 0 0.50 NA NA 0 0.0769231 0
Hydrophis schistosus 0 0 0.02 0.05 0.0555556 0 0.0000000 0

1.3 Sex ratios

# Plotting sex ratios across years

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  group_by( Species, Year, Month)%>%
  summarise(N = n(),
            females = sum(Sex == "Female"),
            prop.female = females/N)%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month)%>%
  ggplot(aes(Month, prop.female))+
  geom_col(width = 0.5, col = "black")+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  facet_grid(Year ~ Species)

Proportion of females to entire population stays constant over the sampling period.

snakes%>%
  group_by(Species)%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  summarise(females = sum(Sex == "Female"),
            N = n())
Species females N
Hydrophis curtus 102 212
Hydrophis schistosus 304 618
# is proportion of females different from 0.5?
snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  group_by(Species)%>%
  summarise(females = sum(Sex == "Female"),
            N = n())%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$females, .$N, p = 0.5)),
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()
Species estimate statistic p.value parameter conf.low conf.high method alternative
Hydrophis curtus 0.4811321 0.2311321 0.6306857 1 0.4125066 0.5504525 1-sample proportions test with continuity correction two.sided
Hydrophis schistosus 0.4919094 0.1310680 0.7173273 1 0.4518628 0.5320580 1-sample proportions test with continuity correction two.sided

Sex ratio is not different from 0.5 p = 0.71

2 Observed reproductive cycle of H. schistosus

# percentage of gravid females in sample

HS%>%
  count(Gravid)%>%
  mutate(N = sum(n))%>%
  mutate(prop.gravid = n/N)%>%
  filter(Gravid == "Yes")
Gravid n N prop.gravid
Yes 88 967 0.0910031

Ì¥

# checking the number of gravid females per year

HS%>%
  group_by(Year)%>%
  filter(Gravid == "Yes")%>%
  count(Gravid)%>%
  spread(Gravid, n)
Year Yes
2016 1
2018 75
2019 12

Proper sampling was only conducted in 2018/19 and hence only this period is used for analysis of reproductive cycles.

# calculating the proportion of gravid females per month

gravid <- HS%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  filter(Year == "2018" | Year == "2019")%>% # only for 2018/19
  group_by(Month)%>%
  summarise(N = n(),
            gravid = sum(Gravid == "Yes"),
            prop.gravid = gravid/N)

print(gravid)
## # A tibble: 8 x 4
##   Month        N gravid prop.gravid
##   <fct>    <int>  <int>       <dbl>
## 1 January    102     12      0.118 
## 2 February   134     34      0.254 
## 3 March      129     24      0.186 
## 4 April       54      7      0.130 
## 5 May         36      1      0.0278
## 6 October     16      0      0     
## 7 November    69      3      0.0435
## 8 December    34      6      0.176

Describing change in the relative proportions of gravid females across months and years of sampling.

# plotting prop gravid per month

gravid%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month, fill = list(prop.gravid = 0))%>%
  ggplot(aes(Month, prop.gravid))+
  geom_point(size = 3)+
  geom_path(aes(group = 1), size = 1, linetype = "dashed")+
  geom_segment(aes(x = "April", xend = "April", y = 0 , yend = 0.02), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(aes(x = "April", y = 0.04, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 0.20, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Proportion of gravid females")

Pregnancy from Novemeber to May. Peak in Feb. Observed two live births in April.

  • Anova/GLM with change in proportion of gravid females over time.

3 Development of embryos and eggs

3.1 Sample size

embryos%>%
  summarise(n.mothers = length(unique(Field.Code)),
            n.embryos = length(unique(Embryo.Code)))
n.mothers n.embryos
29 235

3.2 Summary

embryos%>%
  select(Egg.Length..mm.., Egg.Width..mm.., Egg.Weigth..g.., Snout.to.Vent..cm., Embryo.Weight..g.)%>%
  skimr::skim()%>%
  skimr::yank("numeric")%>%
  select(-hist)

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100
Egg.Length..mm.. 80 0.66 39.71 11.72 2.23 31.66 37.82 47.54 73.98
Egg.Width..mm.. 80 0.66 32.23 9.01 2.02 25.56 33.72 37.64 53.09
Egg.Weigth..g.. 51 0.78 14.30 5.45 4.23 10.75 14.00 17.25 31.00
Snout.to.Vent..cm. 111 0.53 17.72 5.24 6.30 14.38 16.40 21.08 30.40
Embryo.Weight..g. 66 0.72 4.83 6.28 0.00 0.00 3.00 8.00 20.44
embryos%>%
  group_by(Field.Code)%>%
  filter(Sex != "")%>%
  summarise(prop.female = sum(Sex == "Female")/n())%>%
  skimr::skim(prop.female)%>%
  skimr::yank("numeric")%>%
  select(-hist)

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100
prop.female 0 1 0.53 0.31 0 0.31 0.46 0.71 1
embryos%>%
  filter(Sex != "")%>%
  summarise(females = sum(Sex == "Female"),
            N = n())
females N
86 166
broom::tidy(prop.test(86, 166, p = 0.5))%>%
  select(estimate:conf.high)
estimate statistic p.value parameter conf.low conf.high
0.5180723 0.1506024 0.6979603 1 0.4395567 0.5957383

The sex ratio in clutches is not significantly different from 0.5.

3.3 Infertility

embryos%>%
  count(Embryo)%>%
  mutate(N = sum(n), rate = n*100/N)
Embryo n N rate
Absent 5 235 2.12766
Present 230 235 97.87234

3.4 Matrotrophy

Do female H. schistosus expend energy in the development of embryos? or Are the eggs formed and yolk only provides nourishment?

embryos%>%
  ggplot(aes(Embryo.Weight..g., Yolk.Weight..g.))+
  geom_point(size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)

embryos%>%
  select(Yolk.Weight..g., Embryo.Weight..g.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Yolk.Weight..g. ~ Embryo.Weight..g., data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 13.0905565 0.2965229 44.14687 0 0.4071222
Embryo.Weight..g. -0.4526283 0.0450509 -10.04705 0 0.4071222

Yolk weight decreases as embryo weight increases.

embryos%>%
  ggplot(aes(Embryo.Weight..g., Egg.Weigth..g..))+
  geom_point(aes(col = Yolk.Weight..g.), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_color_viridis(name = "Yolk Weight (g)")+
  labs(x = "Embryo Weight (g)", y = "Total Egg Weight (g)")

embryos%>%
  select(Egg.Weigth..g.., Embryo.Weight..g.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Egg.Weigth..g.. ~ Embryo.Weight..g., data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 13.5441434 0.4114909 32.91481 0 0.2188866
Embryo.Weight..g. 0.3676922 0.0557912 6.59050 0 0.2188866

Egg weight increases as embryo weight increases

Decrease of fat bodies as embryos develop indicates matrotrophic nutrition.

  • Photos of females dissected ordered according to embryo dev stage

4 Reproductive effort

# Cleaning reproductive effort data

re <- embryos%>%
  left_join(snakes, by = c("Date", "Field.Code", "Species"))%>%
  select(Date, Field.Code, Embryo.Code, Species, Egg.Length..mm..:Tail.Length..cm..x,
         Egg.Weigth..g..:Sex.x, Snout.to.Vent..cm..y:Tail.Length..cm..y, Weight..g.,
         -Body.Length..cm..x, - Body.Length..cm..y)%>%
  filter(Species == "Hydrophis schistosus")%>%
  rename(Embryo.SVL = Snout.to.Vent..cm..x,
         Embryo.TL = Tail.Length..cm..x,
         Embryo.Sex = Sex.x,
         Female.SVL = Snout.to.Vent..cm..y,
         Female.TL = Tail.Length..cm..y,
         Mother.Wt = Weight..g.
         )%>%
  filter(Female.SVL > 50)

4.1 Increase in clutch size with female age

Does the amount of enery expended by female H. schistosus reduce with age?

re_clutch <- re%>%
  select(Field.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL)%>%
  group_by(Field.Code)%>%
  summarise(clutch.size = n(),
            Clutch.wt = sum(Egg.Weigth..g..),
            Mother.wt = last(Mother.Wt),
            Female.SVL = last(Female.SVL))%>%
  mutate(rcm = Clutch.wt/(Mother.wt - Clutch.wt))
re_clutch%>%
  ggplot(aes(Female.SVL, clutch.size))+
  geom_point(size =3)+
  geom_smooth(method = "lm")

re_clutch%>%
  select(Female.SVL, clutch.size)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(clutch.size ~ Female.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) -12.0816516 4.4876375 -2.692208 0.0154238 0.5466726
Female.SVL 0.2134755 0.0471483 4.527747 0.0002976 0.5466726

The number of eggs borne by females increases with age (SVL).

4.2 Change in overall reproductive effort with age

re_clutch%>%
  ggplot(aes(Female.SVL, rcm))+
  geom_point(aes(col = clutch.size), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed")+
  scale_y_continuous(name = "Relative clutch mass")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Clutch size")+
  theme(legend.text = element_text(size = 10))

re_clutch%>%
  select(Female.SVL, rcm)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(rcm ~ Female.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 0.6886163 0.4077131 1.6889727 0.1170168 0.063695
Female.SVL -0.0038362 0.0042458 -0.9035133 0.3840327 0.063695

The overall reproductive effort does not change with female age.

4.3 Change in reproductive effort per embryo with female age

Does the effort per embryo change with female age?

# reproductive effort per embruo

re_embryo <- re%>%
  select(Field.Code, Embryo.Code, Egg.Weigth..g..,  Mother.Wt, Female.SVL, Embryo.Sex, Embryo.SVL)%>%
  group_by(Field.Code)%>%
  mutate(clutch.wt = sum(Egg.Weigth..g..), 
         Female.wt = Mother.Wt - clutch.wt)%>%
  group_by(Field.Code, Embryo.Code)%>%
  summarise(inv = Egg.Weigth..g../Female.wt, 
            Female.SVL = last(Female.SVL),
            Embryo.Sex = Embryo.Sex,
            Embryo.SVL = Embryo.SVL)%>%
  ungroup()
# using residuals to control for embryo development

emsvlinv <- lm(inv ~ Embryo.SVL, data = re_embryo)

re_embryo%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")

re_embryo%>%
  select(Female.SVL, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(inv ~ Female.SVL + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:p.value, r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 0.0858802 0.0260524 3.296444 0.0017188 0.4433341
Female.SVL -0.0007894 0.0002588 -3.050593 0.0035109 0.4433341
Embryo.SVL 0.0023248 0.0004171 5.573820 0.0000008 0.4433341

The relative egg mass (controlled for embryo development) reduces with female SVL.

4.4 Difference in reproductive effort for male and female embryos

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid, shape = Embryo.Sex))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  select(Embryo.Sex, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(inv ~ Embryo.Sex + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:p.value, r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 0.0165286 0.0090610 1.824151 0.0735636 0.3859103
Embryo.SexMale -0.0083114 0.0045803 -1.814613 0.0750384 0.3859103
Embryo.SVL 0.0022764 0.0004455 5.109887 0.0000042 0.3859103

5 Mortality in bycatch

5.1 Overall mortality

snakes%>%
  group_by(Species)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species n N prop.dead
Hydrophis curtus 120 282 0.4255319
Hydrophis schistosus 174 967 0.1799380
snakes%>%
  group_by(Species)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)%>%
  ungroup()%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)),
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(estimate1:p.value)
estimate1 estimate2 statistic p.value
0.4255319 0.179938 71.81006 0

5.2 Mortality by age class

snakes%>%
  filter(!is.na(Class))%>%
  group_by(Species, Class)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species Class n N prop.dead
Hydrophis curtus Adult 40 62 0.6451613
Hydrophis curtus Juvenile 41 125 0.3280000
Hydrophis schistosus Adult 138 648 0.2129630
Hydrophis schistosus Juvenile 5 105 0.0476190
Hydrophis schistosus Neonate 3 8 0.3750000
snakes%>%
  filter(!is.na(Class))%>%
  group_by(Species, Class)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)),
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(estimate1:p.value)
Species estimate1 estimate2 statistic p.value
Hydrophis curtus 0.6451613 0.328000 15.71186 0.0000738
Hydrophis schistosus 0.2129630 0.047619 17.68174 0.0001447

5.3 Mortality by sex

snakes%>%
  filter(Sex != "")%>%
  group_by(Species, Sex)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species Sex n N prop.dead
Hydrophis curtus Female 36 102 0.3529412
Hydrophis curtus Male 52 110 0.4727273
Hydrophis schistosus Female 62 304 0.2039474
Hydrophis schistosus Male 72 314 0.2292994
snakes%>%
  filter(Sex != "")%>%
  group_by(Species, Sex)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)),
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(estimate1:p.value)
Species estimate1 estimate2 statistic p.value
Hydrophis curtus 0.3529412 0.4727273 2.6538718 0.1032980
Hydrophis schistosus 0.2039474 0.2292994 0.4448479 0.5047918

5.4 Mortality by female reproductive status

HS%>%
  filter(Sex == "Female",
         Class == "Adult")%>%
  group_by(Gravid)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Gravid n N prop.dead
41 151 0.2715232
Yes 16 85 0.1882353
HS%>%
  filter(Sex == "Female",
         Class == "Adult")%>%
  group_by(Gravid)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)%>%
  ungroup()%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)),
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(estimate1:p.value)
estimate1 estimate2 statistic p.value
0.2715232 0.1882353 1.629856 0.2017229